function [ F ] = getDeformationGradientFD( x, y, domain2D, theta )
%GETDEFORMATIONGRADIENTFD Summary of this function goes here
%   Detailed explanation goes here

epsilon = 0.001;

nE = size(domain2D.e1, 1);

F = zeros(2, 2);

% up

up = zeros(2, 1);
p = [x, y+epsilon];

for j = 1 : nE        
        
    qi0 = domain2D.e1(j);
    qi1 = domain2D.e2(j);
        
    q0 = [domain2D.vx(qi0), domain2D.vy(qi0)];
    q1 = [domain2D.vx(qi1), domain2D.vy(qi1)];
              
%     qn = [domain2D.nx(j), domain2D.ny(j)];
        
    v0 = evaluateQuadraturePoint(p, q0);
    v1 = evaluateQuadraturePoint(p, q1);
        
    up(:) = up(:) + ((v0 + v1) * domain2D.l(j) / 2.0 * theta(2*j-1:2*j)); 
end

% down

down = zeros(2, 1);
p = [x, y-epsilon];

for j = 1 : nE        
        
    qi0 = domain2D.e1(j);
    qi1 = domain2D.e2(j);
        
    q0 = [domain2D.vx(qi0), domain2D.vy(qi0)];
    q1 = [domain2D.vx(qi1), domain2D.vy(qi1)];
              
%     qn = [domain2D.nx(j), domain2D.ny(j)];
        
    v0 = evaluateQuadraturePoint(p, q0);
    v1 = evaluateQuadraturePoint(p, q1);
        
    down(:) = down(:) + ((v0 + v1) * domain2D.l(j) / 2.0 * theta(2*j-1:2*j)); 
end

% left

left = zeros(2, 1);
p = [x-epsilon, y];

for j = 1 : nE        
        
    qi0 = domain2D.e1(j);
    qi1 = domain2D.e2(j);
        
    q0 = [domain2D.vx(qi0), domain2D.vy(qi0)];
    q1 = [domain2D.vx(qi1), domain2D.vy(qi1)];
              
%     qn = [domain2D.nx(j), domain2D.ny(j)];
        
    v0 = evaluateQuadraturePoint(p, q0);
    v1 = evaluateQuadraturePoint(p, q1);
        
    left(:) = left(:) + ((v0 + v1) * domain2D.l(j) / 2.0 * theta(2*j-1:2*j)); 
end

% right

right = zeros(2, 1);
p = [x+epsilon, y];

for j = 1 : nE        
        
    qi0 = domain2D.e1(j);
    qi1 = domain2D.e2(j);
        
    q0 = [domain2D.vx(qi0), domain2D.vy(qi0)];
    q1 = [domain2D.vx(qi1), domain2D.vy(qi1)];
              
%     qn = [domain2D.nx(j), domain2D.ny(j)];
        
    v0 = evaluateQuadraturePoint(p, q0);
    v1 = evaluateQuadraturePoint(p, q1);
        
    right(:) = right(:) + ((v0 + v1) * domain2D.l(j) / 2.0 * theta(2*j-1:2*j)); 
end

px = (right - left) / (2 * epsilon);
py = (up - down) / (2 * epsilon);

F(1,1) = 2 * px(1);
F(2,2) = 2 * py(2);
F(1,2) = px(2) + py(1);
F(2,1) = F(1,2);

end

